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A new tool for visualization and analysis of system dynamics is introduced: the 
phasegram. Its application is illustrated with both classical nonlinear systems 
(logistic map and Lorenz system) and with biological voice signals. Phasegrams 
combine the advantages of sliding-window analysis (such as the spectrogram) 
with well-established visualization techniques from the domain of nonlinear 
dynamics. In a phasegram, time is mapped onto the x-axis, and various vibra- 
tory regimes, such as periodic oscillation, subharmonics or chaos, are identified 
within the generated graph by the number and stability of horizontal lines. A 
phasegram can be interpreted as a bifurcation diagram in time. In contrast to 
other analysis techniques, it can be automatically constructed from time- 
series data alone: no additional system parameter needs to be known. Phase- 
grams show great potential for signal classification and can act as the 
quantitative basis for further analysis of oscillating systems in many scientific 
fields, such as physics (particularly acoustics), biology or medicine. 
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1. Introduction 

Oscillating systems, either forced or self-sustaining, are found in many branches 
of physical or biological science. They range from simple harmonic oscillators to 
complex nonlinear systems. Particularly in signals (i.e. time-series data) from 
biological systems, three principal modes of operation are frequently observed: 
periodic oscillation, subharmonics and chaos [1-4]. 1 

In order to assess the complexity of a dynamic system, several analysis para- 
meters have been developed, such as the correlation dimension [5] or the 
Lyapunov exponent [6]. These measures work well on stationary signals 2 
(where system parameters are not varying in time), but are less well equipped 
to deal with signals where system parameters change — as for most (biophysi- 
cal data [7,8]. To cater for these needs, analysis can be performed with sliding 
windows (i.e. the progressional analysis of shorter portions of the signal, 
extracted at consecutive time instants). 

In many applications, the spectrogram is the de facto standard for continuous 
windowed analysis of periodicity and /or the dynamic evolution of (bio)physical 
time-series data. A spectrogram is a type of sliding power spectral analysis. 
When created with a narrow-band setting (resulting in a high-frequency resol- 
ution), the (dis)appearance of spectral peaks in time, owing to bifurcations 
occurring in the underlying system, provides certain important insights into 
the system's dynamics [8,9]. 

While the spectrogram is an analysis tool that is easily accessible to a broad 
range of researchers, it does not provide direct information on the system's 
dynamics in phase space, i.e. a mathematical space where all possible states 
of a system are represented (see Nolte [10] for an essay on the history of 
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Figure 1. Illustration of limitations of both spectrograms and phase space diagrams. Left column: I. deterministic chaos, created by a logistic map; middle column: 
II. stochastic signal created by additive synthesis with randomized phases; right column: III. signal created by a logistic map where the parameter a was gradually 
varied from 3.5375 to 3.6. [a] All 1000 data points of each time series, (b) Amplitude spectrum, (c) Spectrogram, (d) Phase space embedding. Note that both signals 
I. and II. look identical when analysed by Fourier series {b) and (c), whereas their different nature becomes apparent in phase space (cf). The additional complexity of 
signal III. (period doubling cascade to chaos) is not apparent in the phase space diagram (cf). 



phase space). In phase space, each possible state of the system 
corresponds to one unique point. Temporal developments of 
the system are delineated by the so-called trajectories within 
the phase space. For mechanical systems, each degree of 
freedom (d.f.) of the system is mapped onto an individual 
dimension in n-dimensional phase space. 

In the following paragraphs, we will show the advantages 
and disadvantages of the discussed visualization methods, 
i.e. sliding-window analysis (represented by the spectrogram) 
and phase space diagrams. 

The spectrogram's inability to distinguish deterministic chaos 
from a random process is illustrated in figure 1. A synthesized 
signal was generated from a logistic map with parameter a set 
to a value of 3.6 (please refer to case 1, later in the text, for a defi- 
nition of the logistic map) at a sampling rate of 1000 Hz. The 
resulting irregular time-series data are shown in the left column 
of figure la (signal I). The signal's spectrum is displayed in 
figure lb, and a spectrogram is shown in figure lc. The phase 
space diagram 3 (figure Id) reveals that the system's trajectory 
is aligned along a parabola [4, p. 357]. 

The middle panel (II) of figure 1 shows raw and analysis 
data of a series of superimposed sine waves with randomized 
phase offsets, whose amplitudes are derived from the spec- 
tral analysis of the logistic map signal (I). Although signals 
I and II are virtually indistinguishable by spectral analysis 
(figure lb,c, left and middle column), the phase space for 



signal II contains no pronounced attractor but a random 
distribution, revealing the stochastic nature of signal II. 

In order to illustrate the limitation of phase space dia- 
grams, a signal exhibiting a series of bifurcations has been 
created from the logistic map, by gradually varying parameter 
a from 3.5375 to 3.6 (signal III in figure 1, right column). The 
development of the period doubling cascade is clearly visible 
in the spectrogram (being a sliding-window analysis that is 
capable of visualizing temporal developments), but the fact 
is obscured in the phase space diagram (figure Id, right 
column), where all of the different regimes are superimposed. 

Non-stationary signals, potentially exhibiting multiple 
bifurcations over the course of a recorded time series, occur 
frequently in the (bio)physical domain, and they are particu- 
larly common in acoustical data. Here, we address the need 
for an intuitive visualization tool that displays the time- 
varying dynamics of these systems. We introduce a method 
that combines the advantages of both sliding-window analy- 
sis (i.e. the sensitivity to temporal developments in a signal) 
and phase space diagrams (i.e. their close relation to the 
system's underlying dynamics). This new tool, termed the 
'phasegram', is able to visualize system dynamics over time 
in a single two-dimensional graph. The usefulness of this 
method will be exemplified for the specific field of voice 
science using a series of examples with different levels of 
complexity and control over the underlying system. 



2. Methods 

2.1. Phasegram generation 

The phasegram generation process is analogous to creating 
electroglottographic (EGG) wavegrams, a method previously 
developed by Herbst et at. [11]. Phasegram generation is outlined 
below (see also figure 2), and will be described in more detail in 
the following paragraphs. 

— Generation of two-dimensional phase portraits, extracted 
from consecutive windows in time (optional: creation of a 
so-called phase portrait movie). 

— Creation of Poincare sections through the two-dimensional 
phase portraits. 

— Conversion of Poincare section crossings into trajectory 
histograms. 

— Conversion of trajectory histograms into 'trajectory strips' by 
colour-coding. 

— Combination of rotated 'trajectory strips' into the final 
phasegram. 

2.2. Step 1: generation of phase portraits, extracted at 
consecutive points in time 

A phase portrait is the geometrical representation of the trajec- 
tories of a dynamical system, providing a description of the 
system's evolution in time [9]. If the governing differential 
equations of a dynamical system are not known (the typical 
case when analysing a biological system), then the dynamics of 
the phase space can be analysed via attractor reconstruction 
[12,13]: the attractor is defined as a set on the phase plane to 
which all neighbouring trajectories converge [4]. In attractor 
reconstruction, a two-dimensional vector 

x(t) = (B(t),B(t + r)), where t > 0, (2.1) 

is defined, based on the analysed signal. The time series B(t) 
appears as a trajectory x(t) in a two-dimensional phase space [4, 
p. 438]; in other words, the signal is projected against a delayed ver- 
sion of itself (see the electronic supplementary material, figure SI). 
This time delay r should typically be in the range of one-tenth to 
one-half the mean orbital period around the attractor [4, p. 440]. 
A suitable value for r is typically found when the analysed signal's 
autocorrelation function first passes through zero. 

An alternative method for determining the ideal value for r 
has been proposed by Frazer & Swinney [14], who suggest con- 
sidering the first minimum in the mutual information function of 
the attractor as the proper time delay r. This approach should 
ensure the minimized redundancy of information between the 
embedding axes. 

As a third option, the analysed time-series data can be converted 
into an analytic signal by means of a Hilbert transformation. The 
Hilbert transform shifts the phase of each negative frequency com- 
ponent of the analysed signal by +90° ( 77/2 radians) and the phase 
of the positive frequency components by -90° (—tt/2 radians). For 
a purely sinusoidal signal that contains only one frequency com- 
ponent, for instance, applying the Hilbert transform for phase 
portrait generation would have the same effect as using a lag of 
one quarter of the sine wave's period. In order to create a phase por- 
trait from Hilbert-transformed data, the real values of the analytic 
signal are plotted against the imaginary values [15]. This approach 
is particularly useful for longer signals with either highly variable 
cycles or without apparent fundamental frequency. 

While attractor reconstruction can theoretically be performed in 
an unlimited number of dimensions, two dimensions are chosen 
for practical reasons for the purpose of phasegram generation. 
In a discrete-time (i.e. sampled) signal, equation (2.1) becomes 

y[z] = (x[z],x[z + n]), n>0, (2.2) 



where i is the sample index and n is the number of samples delay 
between the two versions of the signal. The process of phase por- 
trait generation is further detailed in the electronic supplementary 
material. 

For the creation of a phasegram, a series of phase portraits is 
required. These are constructed for consecutive portions of the 
analysed signal, centred at constantly progressing time instants. 
These portions are rectangular windows of the signal, defined 
as [x[i — m] . . . x[i + m]], where m is half the window length, 
and i is the sample index around which the window is centred. 
For the consecutive generation of multiple phase portraits, i is 
advanced by steps of / s t, where / s is the sampling frequency 
and r is the time step (usually in the range of 0.01-0.02 s). 

For optimal phasegram output, any DC offset should be 
removed from the analysed time-series data before phase portrait 
generation. 4 In extreme cases where longer signals exhibit a 
considerable baseline drift, a phase-preserving high-pass 
filter with a cut-off frequency of 1 or 2 Hz can be applied. 

As an optional visualization tool, a 'phase portrait movie' 
may be created, by converting the phase portraits (extracted at 
time instants increasing by 0.04 s, which corresponds to the 
default video frame rate of 25 fps) to a video file. The original 
signal can be included as an audio channel, by converting the 
signal to an audio file (WAV file type) with a sampling frequency 
of Fs = 1/dt (Hz), where dt (given in seconds) is the time inter- 
val at which the original time-series data were sampled. The 
synchronous playback of both the auditory stimulus and time- 
varying phase portraits provides valuable insights into the 
analysed sequence. It can also serve as the basis upon which 
the manual rotation of Poincare sections (see below) can be per- 
formed, if so desired. The phase portrait movies for the examples 
discussed in this manuscript are available as electronic 
supplementary material. 

2.3. Step 2: creation of Poincare sections through the 
two-dimensional phase portraits 

A Poincare section is the intersection of a dynamic system's tra- 
jectory in the phase space with a certain lower dimensional 
subspace [9, p. 64]. In a system with n-dimensions, the Poincare 
section is an (n — l)-dimensional surface of section [4, p. 278]. 
Hence, the Poincare section of a two-dimensional phase 
portrait, as created by means of attractor reconstruction, is a 
(one-dimensional ) line . 

The Poincare section generation process is illustrated in 
figure 2d: a line is drawn at a certain angle (horizontal in the case 
of figure 2d) through the two-dimensional phase portrait. The 
number and location of intersections of the phase space trajectory 
and the line (the red dots in figure 2d) are determined and stored 
for further analysis. For the purpose of phasegram generation, the 
Poincare section is made through the entire phase portrait (see 
the electronic supplementary material, movie SI and the audio 
track contained therein for an illustration of the principle). 

The default orientation of the Poincare sections through the 
phase portraits is either horizontal or vertical. In some cases, 
these angles cannot capture crucial aspects of the emerging 
attractor. In such cases, the angle of the Poincare sections (i.e. 
the lines drawn across the phase portraits) needs to be adjusted. 
An optimal angle can be determined visually by inspection of the 
phase portrait movies. Alternatively, a fully automated auto- 
nomous algorithm to determine the rotation angle of the 
Poincare section can be used, which is described in the electronic 
supplementary material. This algorithm aims to find the Poincare 
sections that are best suited to reveal the full complexity of 
the analysed signal's attractor. The phase portrait angles of 
the examples shown in the subsequent case studies were all 
determined with the suggested algorithm. 
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Figure 2. Process of phasegram generation, (a) Spectrogram of a signal synthesized from the logistic map (see text), (b) Time-domain signal of x[i]. (c) Three 
portions of the signal, extracted at times t = 1.35, t = 4.9 and t = 8.57 s. (rf) Phase portraits from the above signals, created by attractor reconstruction. 
A Poincare section was created along the x-axis (orange), yielding intersection points with the trajectory (red dots), (e) Histograms of trajectory intersection 
points with Poincare sections for all three extracted signal portions. For better visibility, a very large histogram bin width of 0.025 was chosen, (f) 'Trajectory 
strips': colour-coded histograms of Poincare sections through phase portraits (see text); (g) phasegram from the signal displayed in (a) and (b). The markers 
at t = 1.35, t = 4.9 and t = 8.57 s represent the position of the three trajectory strips from (f) within the graph. 



The effect of the Poincare section angle in phasegram gener- 
ation is exemplified in figure 3 by a signal generated with a 
Lorenz system [16, equations 25-27] — for details see below. 



Poincare sections were created at 0° (left panel) and 45° (right 
panel). The effect of the angle variation is seen in the histograms, 
the trajectory strips and the resulting phasegrams (figure 3f-h). 
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Figure 3. Effect of the Poincare section angle on phasegram generation, illustrated with a Lorenz system, (a) Parameter B of the Lorenz system, varied from 250 to 28 
(see text), (b) Spectrogram of generated time-domain signal (output variable x of the Lorenz system), (c) Time-domain signal, (d) A 200 ms portion of the time-domain 
signal, centred around t = 9 s. (e) Phase portraits from the above signals, created by attractor reconstruction. A Poincare section was created along the x-axis (left panel) 
and at an angle of 0.35 tt radians, as determined by the algorithm described in the electronic supplementary material (right panel), yielding intersection points with the 
trajectory (red dots). ( f ) Histograms of the trajectory intersection points on Poincare sections in the two conditions, (g) 'Trajectory strips': colour-coded histograms with 
Poincare sections through phase portraits, (h) Phasegrams from the signal displayed in panel (c) for both Poincare section angles, respectively. The markers at t = 9 s 
represent the position of the trajectory strips from (g) within the graph. See text for the note on the period doubling sequence at t = 3.4-3.8 s in the right panel. 



Note that the phasegram created with the algorithmic angle 
selection (right panel in figure 3h) reveals a subharmonic 
sequence around t = 3.4-3.8 s, which is not apparent in the 
phasegram created with an arbitrary angle of 0 radians (left 
panel in figure 3h) — see also electronic supplementary material, 
movie S2. 

2.4. Step 3: conversion of Poincare sections into 
trajectory histograms 

For each phase portrait, a histogram of the trajectory intersection 
points is generated (figure 2e). The histogram bin width is depen- 
dent on the overall signal amplitude, and on the graph height of 
the resulting phasegram. For the phasegrams presented in this 
manuscript, values in the range of 0.001-0.01 were used for 
signals that were normalized to [ — 1 ... 1]. 

2.5. Step 4: conversion of trajectory histograms into 
'trajectory strips' by colour-coding 

The maximum number of trajectory intersection points (hist max ) 
within one single histogram bin from the Poincare section histo- 
grams is determined for all analysed phase portraits. All Poincare 
section histograms are then converted (flattened) to 'trajectory 
strips', i.e. they are colour-coded by 

^ 4°4 99 /IE +i )' (2 - 3) 

where col[k] is the colour intensity (ranging from 0 to 1) and h[k] 
is the number of trajectories found in the A:th histogram bin (see 
figure 2/ for examples). 

2.6. Step 5: concatenation of rotated 'trajectory strips' 
into a phasegram 

All resulting colour-coded 'trajectory strips' are rotated anti- 
clockwise by 90°, and concatenated in sequence to form the 
phasegram (see figure 2g, where the process is exemplified 
by three out of 500 total trajectory strips). In the phasegram, 
time is displayed on the x-axis, representing the instants in 
time at which the individual phase portraits have been created 
(algorithm step 1). The y-axis corresponds to the bin index of 
the Poincare section histogram. The colour intensity shows 
the frequency of occurrence of phase portrait trajectory inter- 
sections with the Poincare section, as seen in the respective 
histograms. 

For the purpose of this manuscript, phasegram generation 
was performed using custom software 5 written by author 
C.T.H. in Python using the matplotlib library [17], and the 
phase portrait movies were created from successive still 
images with the free software ffmepg [18]. 

2.7. Analysed scenarios 

In order to illustrate their applications for yielding insights into 
time-varying biological signals, phasegrams have been gener- 
ated for five scenarios of increasing complexity. The first two 
cases are general systems with well-known dynamics. The 
third case is a computational simulation of vocal fold vibration. 
The last two cases were generated from real signals: the first of 
these was created with an excised larynx set-up in which all the 
crucial underlying parameters are experimentally controlled. In 
the other example, a human singing signal, the oscillating 
system is only controlled in a subjective arbitrary manner by 
the singer ('increasing and decreasing the intensity, trying to 
avoid abrupt changes in sound colour'). 



3. Results 

3.1. Case 1: logistic map. Period doubling cascade 

A synthetic signal was created from the logistic map: 

x[i + l]=ax\i](l-x\i]). (3.1) 

Equation (3.1) was evaluated for values of a in the range of 
[3.4 . . . 3.65], and the change of the parameter a was mapped 
linearly onto a time interval of 10 s with a sampling frequency 
of 1000 Hz. The resulting signal was then up-sampled to 
44 100 Hz with the software Praat [19]. This processing step 
introduced additional data points into the time series by 
means of band-limited interpolation using a sine function 
kernel (sinc(x) = sin(ir x)/(tt x); see [20]), thus allowing display 
of the logistic map time series as limit cycles in the phase por- 
traits (see step 1 and figure 2c ,d)- The DC offset was removed 
from the resulting time series data by simple subtraction of 
the mean, and the signal was normalized to [—1 ... 1]. 

The relationship between a portion of the time-domain 
signal, its phase portrait and its visual appearance in the 
phasegram is clearly seen in figure 2c, d and g. 

— For a sine wave, i.e. the simplest form of periodic vibration, 
the phase portrait is a simple limit cycle, having the form of 
a circle or ellipse — left panel in figure 2d. Because the signal 
is perfectly periodic (figure 2c), only two intersection 
points are seen in the phase portrait, even though the Poin- 
care section is intersected by the trajectory 12 times in each 
direction. This is reflected in the histogram (figure 2e), 
where a total of 24 intersections are distributed over only 
two histogram bins. Consequently, only two horizontal 
lines are seen in the phasegram in the case of periodic oscil- 
lation/the limit cycle (figure 2g, t = 0-2.66 s). 

— In the period doubling case (centre column of panels in 
figure 2), the phase portrait trajectory must complete two 
revolutions in the limit cycle before repeating itself. Because 
there are a total of four intersection points along the Poin- 
care section, the intersections (a total of 23) are distributed 
over four histogram bins (middle panel in figure le). Conse- 
quently, four horizontal lines are seen in the phasegram in 
the interval of t = 2.6-6.23 s in figure 2g. 

— The third column of panels in figure 2 exemplifies the case 
of a non-periodic signal: a complex pattern is seen in the 
phase portrait. A total of 24 intersection points along the 
Poincare section are distributed over a large number (20) 
of histogram bins. Because the signal is non-periodic at 
t = 6.9-9.1 s and t = 9.4-10 s, the histogram bin maxima 
vary greatly over time, and the phasegram obtains a 
'rugged' appearance at these time intervals in figure 2g. 

The resulting phasegram in figure 2g has a striking similarity 
to the standard bifurcation diagram of the logistic map. The 
difference is that the phasegram is derived from real data, with- 
out knowledge of the underlying time-varying parameter. 

3.2. Case 2: Lorenz attractor. Chaos 

The Lorenz system [16, equations 25-27] is defined as: 
X=-aX + oY, 1 
Y = -XZ + rX - Y > (3.2) 
and Z = XY- bZ. 





Figure 4. (a) Schematic of the two-mass model of the vocal folds. Each vocal fold is approximated by two coupled oscillators. The springs and dampers represent 
the viscoelastic properties of the vocal folds. The time-varying lateral displacements xM and xlr were taken as the input to phasegram analysis, (b) Excised larynx 
set-up: sika deer larynx mounted on an air supply tube, EGG electrodes attached on thyroid cartilage at level of vocal folds. 



This system was transformed into a set of difference 
equations: 

+ = <%[;]-*[/]), 1 

y[i + l]=x[t\(r-z[t[)-y[t\ \ (3.3) 
and z[i + 1] = x[i]y[i] — bz[i], J 

and evaluated for a duration of 10 s. The simulation was 
run with a time step of 1/882 s, and a sped-up version of the 
resulting time-series data was then saved with a sampling rate 
of 44 100 Hz. The parameter r was varied gradually from 250 to 
28 over an interval of 7 s, and then kept stable at a value of 28 
for another 3 s (figure 3a). The initial values were set to x = 0, 
y = 20, z = 25. The vector x[i], i = [0 . . . 441 000] was further 
analysed. 

When simulating the Lorenz system (see equations (3.2) and 
(3.3)), parameter a is usually fixed at 10, and b is defined to be 
8/3 [16]. Setting r to a value of 28 results in the well-known 
'strange attractor' of the Lorenz system. A higher value of r 
(above 30) can lead to both stable oscillatory regimes and 
some areas of period doubling bifurcations. This is reflected 
in the phasegram in the right panel of figure 3h: periodic and 
stable from t ~ 0 to 2 s; period doubling from t « 3.5 to 4 s; 
chaotic from t ~ 4 to 10 s. The chaotic nature of the system 
for a stable value of r (t = 7-10 s) is revealed by the phasegram 
in figure 3h (right panel): the intersection points of the phase 
portraits extracted at various instants within this interval vary 
unpredictably, but stay within the regions defined by the 
strange attractor emerging in the phase portrait (see figure 3e). 



3.3. Case 3: self-oscillating two-mass model of the 
vocal folds 

In most mammals (including humans), voice signals are gener- 
ated by the flow-induced, self-sustaining vibration of laryngeal 
tissue [21]. The steady air stream coming from the lungs is 
converted into a time-varying airflow by the oscillation of 
the laryngeal tissue (mainly the vocal folds). The pressure 
variations caused by the time-varying airflow are then propa- 
gated through, and acoustically filtered by, the vocal tract. 
Finally, the result of this process is radiated from the 
mouth (and /or nose) as an acoustical signal [22]. Voice is a 
widely researched physical system that can exhibit a great 
variety of oscillatory behaviour, such as periodic vibration, 



subharmonics, chaos and bifurcations between any of these 
phenomena [3,23-27]. 

In previous research, a simplified two-mass model of voice 
production was created in order to study the effect of asymme- 
tries on vocal fold vibration [28]. In this model, each vocal fold 
is represented by two coupled oscillators (defined by their 
mass, stiffness and damping coefficients; figure 4a). This 
model provides 2 d.f. per vocal fold. It allows for the two 
masses of each vocal fold to oscillate with a phase difference 
(the lower mass typically leading the vibration), thus capturing 
the most essential mechanism of self-sustaining vocal fold 
vibration: the transfer of aerodynamic energy into tissue 
vibrations [29]. The model has the option of simulating the 
effect of asymmetrical vocal fold anatomy/configuration, 
which is a well-known cause of pathological voice production 
and voice nonlinearities [30,31]. 

The two-mass model was run for 10 s with the standard 
parameters used in a previous publication [28], an asymmetry 
coefficient of 0.51, and a time-varying simulated air pressure of 
the lungs ('pressure sweep 7 ) ranging from 0 to 2.5 kPa. The 
time step was 0.05 ms, resulting in a virtual sampling fre- 
quency of 20 kHz. For data analysis, the positions of the 
lower, larger masses were considered. 

Following Steinecke & Herzel [28], the time-varying displa- 
cements of the model's two lower masses were plotted against 
each other (figure 5). Four distinct vibratory regimes were seen 
in the right mass and in the glottal flow (not shown here) for 
the chosen model parameters: periodic, period doubling and 
other subharmonic regimes (see also electronic supplementary 
material, movie S3). Owing to the user-defined asymmetry in 
the model geometry and mechanical properties, the left mass 
had a more complex vibratory pattern than the right 
mass. The right and the left masses are synchronized by 
ratios of 1:1, 5:3, 5:7 and 4 : 2, respectively, for the four 
examples shown in figure 5b-d. 

The four vibratory regimes can be clearly distinguished 
from each other in the phasegram (figure 5e and table 1). 

When relating these four vibratory regimes to the time- 
varying simulated lung pressure used in the simulation 
(figure 5/), a typical hysteresis effect is seen (figure 5g): the 
model's behaviour was different for the pressure increase 
(t = 0-4 s) versus the decreasing pressure (£ = 6-10s). 
The system tended to stay in its current vibratory regime, 
and for certain pressure values (1.6-1.68 kPa; 1.71-1.73 
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Figure 5. Simulation of vocal fold vibration with a simplified two-mass model (Steinecke/Herzel) during a pressure sweep (0-2.5-0 kPa; 0 = 0.51). 
(a) Narrow-band spectrogram of the simulated airflow; (b) and (c) displacement of left (b) and right (c) lower mass of the model as a function of 
time, extracted at t=2, t= 2.76, t = 3.1 and t = 5 s. (d) Phase portraits from the above signals, created by plotting the time-varying position of 
the left lower vocal fold mass against that of the right lower mass. A Poincare section was created at an angle of 0.375 tt radians, yielding intersection 
points with the trajectory (red dots), (e) Phasegram: the vertical markers at f = 0.67; 2.7; 2.91; 3.51; 6.72; 7.27; 7.6 and 9.7 s represent bifurcations, 
i.e. changes from one distinct oscillatory regime to another (see text), (f) Trace of the time-varying simulated lung pressure used to drive the model. 
The vertical markers indicate bifurcation events (see above and text), (g) Bifurcation diagram, showing the distinct vibratory patterns in relation to lung pressure. 
Note the hysteresis caused by the direction of lung pressure change (increasing versus decreasing — see text for details). 



and 1.82-1.88 kPa) more than one vibratory regime was 
possible — see references [31-37] for more detail on hyster- 
esis phenomena in voice. The phonation threshold 



pressure [38], i.e. the minimum pressure required for the 
vocal folds to start and stop the vibration (1.17 and 
1.08 kPa, respectively), was higher for phonation onset 



Table 1. Oscillatory regimes seen in the time series data from the self-oscillating two-mass model simulation of the vocal folds. 



oscillatory regime 


phasegram features 


time (s) 


static 


no oscillation — no line seen in the phasegram (caused by a slight DC offset in the 
analysed signals) 


0 . . . 0.67; 9.7 .. . 10 


periodic 


gradual amplitude variation — two horizontal lines seen in the phasegram 


0.67 . . . 2.7; 7.6 .. . 9.7 


5 '. 3 synchronization 


complex vibratory pattern — six horizontal lines seen in phasegram 


11 1 01 • 7 11 1 & 

LI . . . l.y \, 1 .LI . . . 1.0 


7 : 5 synchronization 


complex vibratory pattern — eight horizontal lines seen in phasegram 


2.91 ...3.51; 6.72... 7.27 


4 : 2 synchronization 


complex vibratory pattern — two horizontal lines in the upper half of the phasegram 
(stemming from the period doubling in the right vocal fold mass — figure 5c) and 
four horizontal lines in the lower half (caused by the period-quadrupling seen in the 
left vocal fold mass — figure Sb) 


3.51 . . . 6.72 



than for phonation offset (representing a subcritical Hopf 
bifurcation). 

3.4. Case 4: excised larynx experiment 

Excised larynx experiments allow experimental investigation 
of vocal production under controlled conditions [39]. The 
larynx (harvested from a freshly deceased individual) is 
mounted on a vertical air tube, and the vocal folds are 
adducted (approximated to the sagittal midline). Humidified 
heated air is blown through the larynx, and the vocal folds 
exhibit flow-induced self-sustaining oscillation if boundary 
conditions are properly set. In such a set-up [40, ch. 1], one 
or more of three basic parameters are usually controlled 
and varied: air pressure, vocal fold adduction and longitudinal 
stress in the vocal folds (vocal fold elongation). 

For the purpose of this study, an excised larynx of a 6.5- 
year-old male sika deer (Cervus nippon) was examined 
(figure 4b). The experimental set-up has been described in 
detail elsewhere [41]. Subglottal air pressure was varied in a 
range of 0-4.1 kPa, as measured with a Keller PR-41X pressure 
sensor positioned 32 cm upstream from the vocal folds. 
Pressure data were captured with a Labjack U6 data acqui- 
sition card at a sampling rate of 1 kHz. The acoustic signal 
was recorded using a DPA 4061 omnidirectional microphone 
positioned 7 cm from the vocal folds. 

Vocal fold vibration was monitored by means of EGG, a 
non-invasive technique that records a physiological correlate 
of vocal fold vibration during phonation [42-44] — see elec- 
tronic supplementary material for details concerning the 
method and the experimental set-up. 

Analysis of the EGG signal from excised sika deer larynx 
oscillations (figure 6a Jo) revealed five distinct vibratory 
regimes: static (no oscillation; not displayed in figure 6b), 
periodic (small amplitude, gradual amplitude variation), irre- 
gular (complex non-periodic signal), period doubling and 
again periodic (stable at larger amplitude). These vibratory 
regimes can be readily distinguished from each other in the 
phasegram (table 2 and figure 6e). 

When relating these vibratory regimes to the time-varying 
simulated lung pressure used in the experiment (figure 6e), a 
hysteresis effect is seen (figure 6/ ), just as in the two-mass 
model (recall figure 5): the system's behaviour was different 
during pressure increase (t = 0-9.5 s) when compared with 
during pressure decrease (t = 9.5-16.99 s). For certain 



pressure values (0.93-1.01; 1.57-1.70; 2.33-2.75 kPa), more 
than one vibratory regime was possible. Contrary to the 
two-mass model example, the phonation threshold pressure 
was lower for phonation onset (0.71 kPa) than for phonation 
offset (0.79 kPa), which is a surprising finding. 

3.5. Case 5: human voice signal 

As for most mammals, human vocalization is created by 
flow-induced self-sustaining oscillation of the vocal folds. In 
non-pathological phonation, four basic vibratory modes 
(called Vocal registers') are observed [21,45]. Of particular 
importance for both speech and singing are two of these: 
the so-called chest register and the falsetto register. In chest 
register, the vocal folds are thickened by the activity of the 
thyroarytenoid muscle (situated within the vocal folds [46]), 
thus introducing a vertical phase delay in the lower versus 
the upper portion of the vocal folds. In general, the duration 
of vocal fold contact (i.e. the (partial) cessation of air flow) 
within one oscillatory cycle is longer in chest than in falsetto, 
resulting in stronger high-frequency harmonics (integer mul- 
tiples of the fundamental frequency) in the generated acoustic 
output, and thus a 'brassier 7 sound [47,48]. In falsetto, the 
thyroarytenoid muscle is relaxed, vertical phase differences 
are less pronounced, the duration of vocal fold contact is 
shorter and weaker high-frequency partials are generated, 
leading to a 'purer', more 'flutey' sound. In speech, both 
males and females mainly use the chest register, which is 
generally lower in fundamental frequency compared with 
falsetto register. In soft- or high-pitched singing (at a 
higher fundamental frequency), the falsetto register is used, 
particularly by females [49]. 

A 52-year-old semi-professional singer (baritone) sang a 
sustained note at vowel /a/ near the upper end of his fre- 
quency range (pitch of C#4, fundamental frequency ca 
277 Hz). He was asked to vary vocal intensity from a mini- 
mum to a maximum, and back to a minimum, without 
taking a breath during the exercise. The goal was to perform 
this manoeuvre on a single musical pitch (i.e. by minimizing 
changes of fundamental frequency) without introducing 
abrupt audible changes into the 'sound colour' of the gener- 
ated sound. Please refer to the electronic supplementary 
material for details on the experimental set-up. 

Human voice production is governed by complex control 
parameters, over which the singer has only partial and 
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Figure 6. Artificial phonation of an excised sika deer larynx produced during a pressure sweep (0-4.1-0 kPa). (a) Narrow-band spectrogram of the electroglotto- 
graphic (EGG) signal; (b) EGG signal, extracted at t = 0.725, t = 2.3, t = 4 and t = 10 s. (c) Phase portraits from the above signals, created by plotting the real 
portion of the Hilbert-transformed EGG signal against its imaginary counterpart. A Poincare section was created at an angle of 0.575 tt radians, yielding intersection 
points with the trajectory (red dots), (d) Phasegram: the vertical markers at t = 0.48, t = 0.98, f = 3.18, t = 6.20, t = 13.14, t = 14.24, t = 15.59 and 
t = 16.07 s represent bifurcations, i.e. changes from distinct oscillatory regime to another (see text), (e) Trace of the time-varying air pressure: the vertical markers 
indicate bifurcation events (see above and text), (f) Bifurcation diagram, showing the distinct vibratory patterns in relation to air pressure. Note the hysteresis 
caused by the direction of air pressure change (increasing versus decreasing — see text for details). 



Table 2. Oscillatory regimes observed in the time series data from a excised sika deer larynx. 



oscillatory regime 


phasegram features 


time 


(s) 




static 


no oscillation — one horizontal line seen in the phasegram 


0... 


0.49; 16.08 . . . 


16.99 


periodic I 


small amplitude, gradual amplitude variation — two horizontal lines seen in the phasegram 


0.49. 


. . 1.00; 15.58 


. . 16.08 


irregular 


non-periodic signal — no stable horizontal lines seen in phasegram 


1.00. 


. . 3.19; 14.30 


. . 15.58 


period doubling 


four horizontal lines seen in phasegram 


3.19. 


. . 6.3; 13.16 . 


. . 14.30 


periodic II 


stable large amplitude — two horizontal lines seen in phasegram 


6.3. 


. 13.16 






(c) 



0.10 
0.05 
0.00 
-0.05 
-0.10 



'falsetto register' 




1.23 



1.23 1.24 
time (s) 



1.24 



time (s) 

'falsetto register' period doubling 
0.2 
0.1 
0.0 
-0.1 
-0.2 
2.14 




2.15 



2.15 2.16 
time (s) 



2.16 



(d) 



0.14 



-0.14 









0.23 


























0 

-0.23 










u 










1 


i i i 


! 


i 


i i 


i 



1.0 
0.5 
0.0 
-0.5 
-1.0 



1.17 



'chest register' 




7.00 



7.01 7.01 
time (s) 



7.02 



-1.17 




-0.14 0 0.14 

Re {x[i\} 



-0.23 0 0.23 

Re {x[i]} 



-1.17 0 1.17 

Re {x[i]} 



(e) 



1.0 



0.5 



-0.5 



-1.0 



1 ! I 
i i 




l ! ! 
i 




■ : ■ 

1 ; 1 




i 
i 








- » r-" - 

i i 




r 
i 






i i 

i : i a 






■ : ■ X 

i i 





0 



10 



time (s) 



abrupt register transition (bifurcation): 
'falsetto' to 'chest' register 

Figure 7. Electroglottographic recording of a 52-year-old singer producing a sustained note (vowel /a/) with intensity variation (soft -loud -soft), (a) Singer's 
attempted intensity of voice production in arbitrary numbers on a nonlinear scale (0, lowest intensity; 1, highest intensity), (b) Narrow-band spectrogram of 
the electrographic (EGG) signal, (c) EGG signal, extracted at t = 1.22, t = 2.14 and t = 7 s. (d) Phase portraits from the above signals, created by plotting 
the real portion of the Hilbert-transformed EGG signal against its imaginary counterpart. A Poincare section was created at an angle of 0.3 tt radians, yielding 
intersection points with the trajectory (red dots), (e) Phasegram: the vertical markers at t = 1.22, t = 2.14 and t = 7 s indicate the time instants at which the 
signals shown in panel (c) were extracted. 



intuitive control. In the example shown in figure 7, the 
singer's intended intensity of voice production was varied, 
attempting to keep all other parameters stable. The plotted 
intensity is a dimensionless quality, expressed on an arbitrary 
nonlinear scale (0, lowest intensity; 1, highest intensity; 
figure 7a). The spectrogram in figure 7b reveals several 
abrupt transitions, suggesting spontaneous changes in the 



underlying voice production mechanism, not intended by 
the singer. They represent unwanted, spontaneous system- 
level behaviour and violate the traditional aesthetic boundary 
conditions of classical singing. The findings, corroborated by 
inspection of the time-domain signal (provided as the audio 
track in electronic supplementary material, movie S5), are 
described in table 3. 



Table 3. Observed system state transitions in the case of an inadequately 
executed intensity variation manoeuvre in singing. 



offset (s) 


observed phenomenon 


0 


periodic 


1.75 


period doubling 


4.2 


chaos 


4.4 


periodic 


4.59 


chaos 


4.63 


periodic 


8.82 


period doubling 


9 


periodic 


9.18 


subharmonic regime 


9.3 


periodic 



EGG signals of the most prominent vibratory regimes were 
extracted at time instants t = 1.22, t = 2.14 and t = 7, and are 
shown in figure 7c. Based on the knowledge gained from pre- 
vious research [48,50,51], the sound production mechanisms of 
the three extracted samples can be classified as 'falsetto regis- 
ter', 'falsetto register with period doubling' and 'chest register', 
respectively. Phase portraits of the three extracted signal por- 
tions are shown in figure Id, and the respective positions in 
time of their Poincare sections are marked as vertical lines in 
the resulting phasegram in figure 7e. The abrupt transitions 
between various vocal fold vibratory regimes as observed in 
the spectrogram (figure 7b) are also clearly reflected in the pha- 
segram. The arrow in figure 7e marks the position of an abrupt 
transition from falsetto to chest register, occurring over an 
interval of ca 18 ms (five vibratory cycles). At that instant of 
time, the fundamental frequency dropped abruptly from 320 
to 242 Hz (as seen in the spectrogram, figure 7b). 

4. Discussion 

Phasegrams provide a new tool to visualize and gain insights 
into system dynamics, particularly for nonlinear dynamical 
systems in which the vibratory regime changes over time. 
In a phasegram, the intuitiveness and accessibility of the 
spectrogram is combined with the more detailed information 
on system dynamics provided by phase space analysis. Basi- 
cally, a phasegram is an empirically derived bifurcation 
diagram in time. However, in contrast to the 'traditional' 
bifurcation map approach, underlying system parameters 
do not need to be known for phasegram generation. This is 
particularly useful in situations where only an output 
signal of the (bio)physical system of interest is available. 

The creation of phasegrams was motivated by the goal of 
visualizing a system's dynamics over time in a single two- 
dimensional graph. To accomplish this, considerable data 
reduction is required: the creation of consecutive time-varying 
two-dimensional phase portraits by means of attractor 
reconstruction results in three-dimensional data. A further 
reduction is achieved by performing Poincare sections through 
the generated phase portraits, thus reducing the data to two 
dimensions: in a phasegram, time is displayed on the x-axis, 
and the time-varying intersection points through the 
Poincare sections are displayed on the y-axis (as an additional 



feature, the frequency of occurrence of the intersections is 
reflected as colour intensity along a virtual z-axis). The resulting 
two-dimensional graph provides intuitive insights into the ana- 
lysed system's dynamic behaviour. Phasegram figures are, due 
to their two-dimensional nature, expected to be particularly 
useful in printed publications. 

For an ideal Poincare section orientation (see 2.3), the 
dynamics of the analysed system can be assessed via the 
number and stability of horizontal lines in the phasegram. 

— One single line: static (no oscillation). 

— In locally stable lines: limit cycle, with optional subhar- 
monics if n > 1, i.e. periodic oscillation for two lines, 
period doubling for four lines, subharmonic regime with 
period tripling [52] for six lines, etc. If the analysed 
nearly periodic signal contains substantial spectral 
energy above the fundamental (i.e. harmonics), additional 
trajectory contours are introduced into the phase por- 
traits, their Poincare sections, and consequently into the 
generated phasegrams. In such cases, low-pass filtering 
the analysed signal may be considered. 

— Multiple, locally unstable lines: irregularity caused by either 
changing system parameters; a quasi-periodic signal with 
more than two individual sinusoids whose frequencies are 
not related to each other by integer ratios; or chaos. 

In the range of examples considered here, we have found 
the phasegram to be a useful tool for classification of (bio)- 
physical signals into three main categories: periodic 
oscillation, subharmonic regimes and chaos. This potential 
needs to be verified in further studies and other systems. 

Just as in other analysis and visualization techniques (e.g. 
spectrograms), the appearance of a phasegram (and hence 
the information content provided by the technique) depends 
on the selection of various parameter values — see electronic 
supplementary material. As is the case in any analysis of 
experimental time-series data, noise in the measurement may 
introduce unwanted effects: when analysing signals with a 
low signal-to-noise ratio, the generated phase portraits and 
hence the phasegram itself can adapt a 'noisy' appearance. 

To generate a phasegram, Poincare sections are needed that 
intersect with a significant portion of the trajectories in phase 
space, in order to reveal the governing dynamics of the 
system. Such an approach cannot guarantee transversality of 
all Poincare sections. In the light of this issue, the applicability 
of the phasegram method depends on the global structure of 
the dynamics of the analysed time series. However, for many 
practical purposes any Poincare section — regardless of trans- 
versality — will produce useful insights into system dynamics. 
Oscillatory dynamics, in general, and acoustic and other 
voice-related physiological data, in particular, are especially 
suited for phasegram analysis: These signals repeatedly cross 
zero by construction (e.g. due to pressure fluctuations), and 
the origin of the phase space is used as the anchor point for 
the Poincare sections. In other cases, where these zero-crossings 
are not present, a simple removal of the DC offset or trend 
found in the data may facilitate a phasegram representation. 

In the past, several methods for two-dimensional space - 
time visualization have been introduced, such as, for example, 
phase space diagrams stemming from numerical solutions of 
equation systems [16], phase space attractor reconstruction 
[12,13], recurrence plots [53] or coarse-grained embeddings of 
time series [54-56]. These tools are very helpful in visualizing 



the dynamics of a system over a certain (and typically short) 
period of time, by analysing the underlying time-series data 
(consisting of a certain number of samples) and generating 
one graph for that particular time window. However, such rep- 
resentations fail to appropriately illustrate changes of system 
state (e.g. from periodic oscillation to period doubling or 
chaos). In particular, visualization tools that do not map time 
onto either the ordinate or the abscissa are unable to pinpoint 
within a single graph the exact points in time at which bifur- 
cations occur. This is particularly crucial for signals that are 
characterized by frequent changes of oscillatory regime, as is 
found frequently in time-series stemming from biosystems or 
acoustic data (see figure 7e for an example). 

To overcome this limitation, we developed the phase- 
gram, based on attractor reconstruction and Poincare 
sections, which are already well-established visualization 
techniques. In contrast to conventional sliding-window 
analysis, such as the spectrogram [57,58], the Wigner-Ville 
distribution [59] or the scaleogram [60], the phasegram 
offers direct insights into the dynamics of the analysed 
system, because it is generated from the signal's embedded 
attractor. No advance knowledge of periodicity or fundamen- 
tal frequency is required for phasegram generation, which 
makes the tool superior to previously presented sliding- 
window visualization approaches such as the time-varying 
sequence of local maxima [24] or EGG wavegrams [11]. 

Although the phasegram is not equipped to measure the 
dimensionality of complex attractors, the examples shown in 
this manuscript suggest that the new method is very well 
suited to visualize the temporal evolution of a dynamical sys- 
tem's attractors, from very simple to complex. Two embedding 
dimensions are sufficient to fully describe (nearly) periodic sig- 
nals and to evaluate the signal's divergence from such a 
periodic pattern. In this respect, phasegrams offer an empirical 
basis for the decision of whether the signal's fundamental fre- 
quency can be calculated, or whether (computer-aided) 
periodicity analysis will give spurious results due to the 



occurrence of subharmonics, quasi-periodicity or other forms 
of irregularity. 

In this paper, the application of phasegrams for system 
dynamics visualization has been exemplified for voice pro- 
duction. However, the technique is by no means limited to 
vocal applications. As our initial logistic map and Lorenz 
system examples show, phasegrams can usefully be applied in 
any branch of physics, biology, chemistry, economics or other 
fields where nonlinear systems are analysed and an intuitive 
indication of system dynamics and how they change over time 
is desired. 
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Endnotes 

lr To make the contents of this manuscript accessible to researchers 
from disciplines other than physics or mathematics, we include 
some tutorial-style background information on these signal types 
as electronic supplementary material. 

2 A signal is classified as non-stationary if its probability distributions 
are time-dependent. Non-stationarity is given if a time-series data 
have a mean, variance or covariance that changes over time. Non- 
stationary features manifest themselves as trends, cycles or random 
walks. Most biological signals (but also, e.g. economic time-series 
or meteorological data) can be approximated as stationary on a 
short scale, but are non-stationary when observed over a longer 
time span. Non-stationarity in (bio)physical systems arises when 
system parameters change, e.g. the neurally controlled tension of 
the vocal folds in voice production. 
3 Alternatively termed 'phase space plot'. 
4 The DC offset is the arithmetic mean of a time-series data. 
5 The phasegram software is available online at the first author's 
website at http : / / www. christian-herbst. org / phasegrams / . 
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